Shock waves in capillary collapse of colloids: a model system for two— dimensional 

screened Newtonian gravity 
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Using Brownian dynamics simulations, density functional theory, and analytical perturbation 
theory we study the collapse of a patch of interfacially trapped, micrometer-sized colloidal particles, 
driven by long-ranged capillary attraction. This attraction is formally analogous to two-dimensional 
(2D) screened Newtonian gravity with the capillary length A as the screening length. Whereas the 
limit A ^ oo corresponds to the global collapse of a self-gravitating fluid, for finite A we predict 
theoretically and observe in simulations a ringlike density peak at the outer rim of a disclike patch, 
moving as an inbound shock wave. Possible experimental realizations are discussed. 
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The dynamics of matter under the influence of long- 
ranged attractions is^ studied intensively in several 
branches of physics [1], in particular with respect to 
inherent instabilities. Most prominently, the structure 
formation in the universe is understood as the conse- 
quence of an instability in self-gravitating matter, and 
cosmological theories in conjunction with numerical sim- 
ulations have been successfully applied to explain the 
dynamical formation of clusters, galaxies, or dark mat- 
ter halos on large scales [Ij. More recently, other sys- 
tems with gravitational-like attractions have been inves- 
tigated, including seemingly unrelated phenomena like 
bacterial chemotaxis [sl, 0] or capillary-driven clustering 
in colloids trapped at fluid interfaces [5, 6]. In these sys- 
tems the interaction is effectively cut off beyond a finite 
range, albeit much larger than the mean interparticle sep- 
aration. 

In a self-gravitating fluid any homogeneous mass dis- 
tribution is unstable with respect to small fluctuations on 
sufficiently large scales [7] (Jeans ^ instability). In systems 
with a cut-off gravitational-like attraction, this instabil- 
ity only occurs below a critical temperature 0, HI • As the 
range of the interaction is scaled [6] from infinity down 
to a microscopic length like the size of the particles, the 
dynamical evolution of the instability crosses over from 
gravitational collapse to spinodal decomposition. A stan- 
dard theoretical approach to the gravitational collapse in 
cosmology is the so-called cold collapse approximation 
(see, e.g., Ref. [8[), within which any force other than 
gravity (in particular the thermal pressure of the fluid) is 
neglected altogether. The view on applications to other 
physical situations raises the natural question how the 
phenomenology of this scenario is affected by a nonva- 
nishing thermal pressure and a large but finite range of 
attractions and specifically how the crossover to the spin- 
odal decomposition scenario occurs. 
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FIG. 1. Proposed sketch of dynamical regimes for a circular 
patch of radius L with particles of radius R as function of 
the range A = A/L of the interaction and a rescaled effective 
temperature T(l + A~^) (see Eq. ^ and below) The rescal- 
ing factor for T leads to a horizontal border separating the 
collapse and dilution regimes. Neglecting a possible temper- 
ature dependence of A, isotherms are parallel to the dashed 
black line. The transition region is bounded approximately 
by the line for which linear stability theory for an infinite ho- 
mogeneous distribution P, ^ predicts that the fastest growing 
density mode has a wavenumber 27v/L. 



Colloidal particles with radii in the micrometer range, 
which are trapped at a fluid interface, lend themselves 
to study this issue. Their weight results in a force / 
on each particle perpendicular to the interface which de- 
forms it and gives rise to long-ranged capillary interac- 
tions between the particles. For large colloid center-to- 
center distances the leading interaction term (dominat- 
ing the collective collapse dynamics) is a pair interaction 
with the potential [9] V{d) = [p/{27r-f)] Ko{d/X) with 
the modified Bessel function Kq^ the capillary length A 
(~ 0(mm)), and the interfacial tension 7. For d < X this 
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FIG. 2. Evolution of the radial density profile for T — 3.1 x 10 ^ (for further parameters, see the main text). Panel (a): 
A = 0.25, comparison between 2D-DDFT (colored lines), and BD (symbols). Panel (b): A = 0.075, only BD (DDFT results 
omitted for clarity) . Unlike as in (a) , small transient peaks distinct from the one at the outer rim are observed due to clustering 
near the center of the collapsing disc. 



reduces to 2D Newtonian gravity, VQ^d) ~ Ind, whereas 
for d > A it decreases exponentially. This soft mat- 
ter system is of particular interest because the range 
A = / {gApm) is tunable via the dependence of 7 on 
temperature T and the concentration of surfactants or 
through the mass density-difference Ap^ between the 
two fluid phases forming the interface; g is Earth's grav- 
itational acceleration. 

Main Results - We study the time evolution of the 2D 
particle number density ^(r = {x, y}^ i) of an initially cir- 
cular patch of radius L with particles of radius R on the 
fiat interface with a homogenous density ^o- The range 
is measured by the dimensionless parameter X = X/L and 
we introduce the effective, dimensionless temperature 



(1) 



For A ^ 00, T is the ratio of the thermal to the (mean- 
field) attractive inner energy of the patch because each 
particle interacts with N ^ QoL"^ other particles. For 
A < L, the number reduces to TV ^ ^O'^^ and this ratio 
approximately equals (1 + X~^)T. Figure [1] summarizes 
qualitatively the dynamic phases we have found. If T is 
large enough (above the black line), the capillary attrac- 
tion cannot confine the particles to the circular patch, 
so that it becomes more and more diluted as time pro- 
gresses. If T is small (dashed line) the patch collapses, 
but the way how this proceeds depends on A. In the 
limit A ^ 00 (gravitational collapse), the fastest growing 
modes span the patch and the evolution is dominated by 
the collapse of the structure as a whole. In the opposite 
limit X ^ R (spinodal decomposition), the fastest grow- 
ing modes have a characteristic length scale well below 
L and the evolution is dominated by the coarsening of 
these domains. We have explored the transition region 



between these limits using perturbation theory, Brownian 
dynamics (BD) simulations, and dynamic density func- 
tional theory (DDFT). Accordingly, the most prominent 
feature of this transition is a density peak forming at the 
outer rim of the collapsing circular patch which exhibits 
properties of a shock wave. The latter is defined by the 
crossing of the characteristic curves of a differential equa- 
tion. For Eq. (|2j) below, in the presence of radial sym- 
metry, these characteristics are the trajectories of rings 
of particles. Rings with initially larger radii travel faster 
than those with smaller radii and, upon crossing, form 
a density singularity as T ^ (Fig. Efa)). This feature 
becomes more pronounced with A becoming smaller, be- 
cause the time for building up the shock wave diminishes 
relative to the collapse time of the whole patch (Fig. [3j). 
On the other hand, for smaller A additional small clus- 
ters form in the interior of the patch and the shock wave 
amplitude is reduced (Fig. [2{b)). 

Theory - We consider a 2D fluid within a mean-field 
approximation appropriate for the long-ranged, capillary 
interactions and take into account nonzero temperature 
and short-ranged interactions through a pressure equa- 
tion of state p{g^T) (Sf, containing the effects of short- 
ranged interactions such as hard or soft cores and sub- 
leading terms in the capillary forces. The particle dynam- 
ics is assumed to be in the overdamped regime with F as 
the associated interfacial mobility of the particles and 
hydrodynamic interactions are neglected. The relevant 
time scale of collapse is Jeans' time T = 7/(r/^^o)- We 
note that the time it takes a particle to diffuse a distance 
L by Brownian motion alone is ~ L^/(r/cBT) = T/T. 
Dimensionless variables are introduced as g = g/go, 
P = p/{kBTgo), r = t/L^ and t = i/T together with 
Eq. (jl]). Mass conservation reads |5| 

§ = -v-(eVC/[e]-rvp), (2) 
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with the dimensionless potential of capihary interaction 
U[g{r)] = l/(27r) J dr' g{r')Ko{\r - r'\/X). For ideahzed 
point particles the pressure is = Q] we have also con- 
sidered a fluid of hard discs of radius described by (TT 
Phdig) = q{Qc^ q)/ {Qc - q), where Qc = (2V3^o^^)" 
the dimensionless density of close packing. Equation 
can be viewed as a simple DDFT for this system [5| 
which, given the long range of the attractions, is expected 
to hold at least for scales much larger than R and to de- 
scribe correctly the collective aspects of the dynamics. 
Deviations are likely to occur for smaller scales and are 
mainly attributable to the local density approximation, 
i.e., the crude form of the term oc Vp for the short-ranged 
interactions. Here, improvement could be achieved using 
more sophisticated expressions from DFT. 

The cold collapse (CC) scenario (formally T = in 
Eq. (|2j)) is best studied using Lagrangian coordinates: the 
characteristic curves of Eq. (j2j) are the radial trajectories 
of infinitesimally thin rings of particles, assuming that 
the initial radial symmetry is preserved by the mean- 
field evolution. If the initial radius is ro, the trajectory 
is described by a time dependent mapping r = rLai'^o^t). 
The Jacobian of this mapping provides the density field 
in Lagrangian coordinates: 
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(ro < 1). 



(3) 



The mapping is well defined and invertible as long as 
two rings do not cross. If two rings cross, the Jaco- 
bian vanishes and the density field exhibits a singular- 
ity. For A ^ oo, the CC approximation yields jH, [llj 
VLa = roVT^ (t < 1) and gLa{ro.t) = 1/(1 -t). Thus, 
the initial homogeneity inside the patch is preserved and 
a singularity arises at t = 1 (i.e., at Jeans' time), when 
all the rings reach the center simultaneously. For finite 
A we have applied perturbation theory in terms of 1/A 
leading to [6] 



rLa ^ ' _ , , , 2 - 47E„ier - rg - ln(4A^) 
' " ^ (4A)2 
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This result predicts the formation of an overdensity for 
the outermost ring (ro = 1) which becomes singular 
(1/^La = 0) at a time ^ 1 + ln(1.44A)/(2A)2 when 
it reaches a radius Vg = rLa(l,^s) ^ 1/(4A). (The sin- 
gularity is a consequence of the CC approximation; it 
is actually regularized into a density peak moving like 
a shock wave by the term oc Vp.) Thus, with decreas- 
ing A the collapse singularity occurs later and at larger 
radii. Intuitively, for decreasing A the range of the in- 
teraction is confined to smaller regions around any point 



of the disc. In addition the outer rim of the disc expe- 
riences no balancing pull from the outside, so that there 
the first overdensity builds up, which then attracts more 
and more particles. 

Numerical Methods and Discussion - In order to test 
this theoretical analysis we have performed simulations 
and solved Eq. ([2]) numerically. Simulation parameters 
were chosen to reflect the conditions in an actual ex- 
perimental realization [5] (patch size L = 1.83 mm, 
particle radius R = 10 /im, capillary potential depth 
/V(27r7) = 0.89 /crT, particle mobility F = 37r7?wateri? 
with T^water being the water viscosity at room tempera- 
ture, and the particle hard-core realized by the repulsive 
part of the Lennard-Jones potential [lH). The simula- 
tions were carried out with N = 1804 particles using 
Brownian Dynamics (BD) [gI, |13| and the radial density 
profile was obtained through angular and ensemble (120 
runs) averages. Numerical solutions were obtained ei- 
ther through a numerical integration of Eq. (|2]) with en- 
forced radial symmetry or through a particle-based (La- 
grangian) integration scheme of the full 2D equation (la- 
beled as 2D-DDFT) [sl,^. In this latter case, the den- 
sity field is probed by a discrete number of (fictitious) 
particles which follow the characteristic curves of Eq. (|2]) 
and the density profile is obtained like the profiles from 
BD simulations. 

We find that the value A = 1.5 (i.e., A ^ 2.7 mm, 
which is the capillary length of the air-water interface 
at ambient conditions) is at the limit of applicability of 
the perturbation theory: the DDFT solutions for T = 
confirm the occurrence of the singularity at the predicted 
values rg ~ 0.16, tg ^ 1.1 and that the formation of the 
singularity indeed slows down and occurs at larger values 
of rg with decreasing A (r^ ^ 0.69, tg ~ 1.7 for A = 0.25). 
The DDFT solutions at nonzero T = 0(10"^ . . . 10"^) 
show that, as for T = 0, an overdensity peak forms at the 
rim while traveling inwards with increasing amplitude, 
but the development into a singularity is inhibited. 

The results from the BD simulations confirm this and 
provide further details about the dynamical evolution, 
see Figs. [2] and [31 For A = 0.25 at T = 3.1 x 10"^, in 
contrast to the collapse at larger A, a ring-shaped densi- 
fied zone forms which is composed of a number of smaller 
clusters. The 2D-DDFT solutions agree well with the 
BD simulation data for all times. This indicates that the 
features we are analyzing are not affected by the details 
of the microscopic correlations, the effect of which can 
be taken into account through the macroscopic pressure 
term. Upon decreasing A further, the formation of even 
more individual clusters inside the disc is observed, as a 
prelude to the spinodal decomposition scenario. For the 
particular value A = 0.075 the averaged radial density 
profile still exhibits the collective behavior of the shock 
wave at the outer rim albeit with a smaller amplitude. 
Additionally, transient peak structures for smaller radii 
become also visible and are eventually absorbed by the 
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FIG. 3. Snapshots from BD simulations for temperature T = 3.1 x 10~^ and A = 1.5 (upper row) or A = 0.25 (lower row). 
Particle distributions in each column have the same radial extent. Clusters (i.e., particles with at least 3 neighbors within a 
distance of 3.25i?) are depicted in red. For A = 1.5 the global collapse appears to be faster than the formation of individual 
smaller clusters. For A = 0.25, small clusters predominantly form at the outer rim and collectively move towards the center. 



shock wave traveling inwards. 

The phenomenology of local clustering due to small 
initial fluctuations is also visible in the azimuthal direc- 
tion, both in the BD simulations and in the 2D-DDFT 
solution. In both cases the initial conditions break the 
radial symmetry, which is expected to be recovered by 
an ensemble average. 

Concerning a possible experimental realization, a suit- 
able value of T can be arranged using colloids with R ~ 
10 /im [5[. Observation of the ring-shaped density build- 
up and the ensuing shock wave requires A/L < 0.25, cor- 
responding to L > 10 mm for an air-water interface. A 
likely relevant issue for experiments is the role of hydro- 
dynamic interactions. We have carried out simulations 
incorporating these on the Rotne-Prager level (as formu- 
lated in Ref. |15]) for the specific case of particles with 
contact angle close to zero (i.e., just touching the inter- 
face). Our results indicate that the qualitative features 
of the evolution discussed here are not affected and the 
collapse is simply accelerated. This is in line with the 
experimentally observed enhanced colloid self-diffusion 
due to hydrodynamic interactions in such a system [l6| . 

Summary and Conclusions - We have studied the 
collapse of a homogeneous, circular patch of colloidal 
particles trapped at a fluid interface by means of an- 
alytical perturbation theory, Brownian dynamics simu- 
lations, and dynamical density functional-like theories. 
The capillary attraction is formally analogous to two- 
dimensional gravity with a tunable cutoff length A. We 
find that a finite value of A strongly influences the col- 
lapse features. While for A ^ oo the evolution is dom- 
inated by the global collapse of the patch, a large but 
finite A induces the formation of a ring-like overdensity 
which quickly becomes singular in the limit of a vanishing 
pressure force (i.e., zero-temperature or cold collapse). 
A nonvanishing pressure regularizes the singularity into 



a collapsing shock wave. System parameters can be cho- 
sen such that these spatio-temporal structures can be 
realized in experiments with micrometer-sized colloids. 
Furthermore, this system appears to be ideally suited to 
investigate the transition from Jeans' gravitational in- 
stability (A oo) to a spinodal instability (A ~ colloid 
radius). 
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